%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the basis matrix of a 1-D spline.
% For the node points outside the approximation interval, the function is 
% approximated as a linear function satisfying the zero-order and 
% first-order continuity at the boundary points. 
% Author:   Xing Guo (xingguo@umich.edu)
% Date of this version: 2017/9/18
% Function: 
% [B,x]=ExtendedSplineBasis(breaks,evennum,k,x,order)
% Input variables: 
%   breaks:     vector of break points
%   evennum:    whether evenly distributed break points
%               0: no
%               others: yes, breaks should only include boundary points
%   k:          polynomial order
%   x:          Evaluation points
%   order:      A scaler or a vector for the evaluation order
% Output variables:
%   B:          Basis matrix (cell if order is a vector)
%   x:          Evaluation points
% Special Note: 
% Uses:     splibas from CompEcon Toolbox by Miranda and Fackler
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [B,x]=ExtendedSplineBasis(breaks,evennum,k,x,order)

% breaks          =   SplineDef.breaks;
% evennum         =   SplineDef.evennum;
% k               =   SplineDef.k;
I_left          =   x<breaks(1);
I_right         =   x>breaks(end);

[B,x]           =   splibas(breaks,evennum,k,x,order);

% Revise the basis matrix for the points on the left of the left boundary
if any(I_left)
    B_left      =   splibas(breaks,evennum,k,breaks(1),[0;1]);
    if length(order)==1
        switch order
            case 0
                B(I_left,:)     =   (x(I_left)-breaks(1))*B_left{2}+...
                                    repmat(B_left{1},[sum(I_left),1]);
            case 1
                B(I_left,:)     =   repmat(B_left{2},[sum(I_left),1]);
            otherwise
                B(I_left,:)     =   0;
        end
    else
        for i=1:length(order)            
            switch order(i)
                case 0
                    B{i}(I_left,:)      =   (x(I_left)-breaks(1))*B_left{2}+...
                                            repmat(B_left{1},[sum(I_left),1]);
                case 1
                    B{i}(I_left,:)      =   repmat(B_left{2},[sum(I_left),1]);
                otherwise
                    B{i}(I_left,:)      =   0;
            end
        end
    end

end
% Revise the basis matrix for the points on the right of the right boundary
if any(I_right)
    B_right         =   splibas(breaks,evennum,k,breaks(end),[0;1]);
    if length(order)==1
        switch order
            case 0
                B(I_right,:)	=   (x(I_right)-breaks(end))*B_right{2}+...
                                    repmat(B_right{1},[sum(I_right),1]);
            case 1
                B(I_right,:)    =   repmat(B_right{2},[sum(I_right),1]);
            otherwise
                B(I_right,:)    =   0;
        end
    else
        for i=1:length(order)            
            switch order(i)
                case 0
                    B{i}(I_right,:)     =   (x(I_right)-breaks(end))*B_right{2}+...
                                            repmat(B_right{1},[sum(I_right),1]);
                case 1
                    B{i}(I_right,:)     =   repmat(B_right{2},[sum(I_right),1]);
                otherwise
                    B{i}(I_right,:)     =   0;
            end
        end
    end
end

